//finite-state, continuous-time, tagging model
//Tim Miller May 2011

TOP_OF_MAIN_SECTION
 //ALLOCATE SPACE IN READ-WRITE MEMORY
 arrmblsize=100000000;
 //gradient_structure::set_GRADSTACK_BUFFER_SIZE(300000000);
 //gradient_structure::set_CMPDIF_BUFFER_SIZE(400000000);
 //gradient_structure::set_NUM_DEPENDENT_VARIABLES(20000);
  time(&start); //this is to see how long it takes to run
  cout << endl << "Start time : " << ctime(&start) << endl; 

GLOBALS_SECTION
  #include <admodel.h>
  #include <time.h>
  time_t start,finish;
  long hour,minute,second;
  double elapsed_time;

DATA_SECTION
 //NO SPACES in (,)
 int ipredict
 int jpredict
 int kpredict
 //number of calendar periods where fishing mortality etc. can change
 init_int n_calendar_periods
 //number of periods of different fishing mortality before it is equal to other fish
 init_int n_nonmixing_periods
 //number of periods of different tag-induced mortality before it is 0
 init_int n_tag_mort_periods
 //number of non-100% reward types
 init_int n_reward_types
 
 //number of regions
 init_int n_R
  
 //the number of release groups for unrecovered fish (could be the number of unrecovered fish...)
 init_int n_release_groups_unrecovered
 //the total number of fish released per release group by state (columns)
 init_imatrix unrecovered_n_released_by_state(1,n_release_groups_unrecovered,1,2*n_R)
 //the number of time intervals with homogenous rates for each release group
 init_ivector n_intervals_per_group(1,n_release_groups_unrecovered)
 //the total number of time intervals with homogenous rates for all release groups
 init_int unrecovered_total_number_of_intervals
 //the duration of each release group time interval
 init_vector unrecovered_intervals(1,unrecovered_total_number_of_intervals)
 //the reward type for each release group
 init_ivector unrecovered_reward_type(1,n_release_groups_unrecovered)
 //the calendar period for each release group time interval
 init_ivector unrecovered_calendar_period(1,unrecovered_total_number_of_intervals)
 //the post-tagging nonmixing period for each release group time interval
 init_ivector unrecovered_nonmixing_period(1,unrecovered_total_number_of_intervals)
 //the post-tagging tag-induced mortality period for each release group time interval
 init_ivector unrecovered_tag_mort_period(1,unrecovered_total_number_of_intervals)
 
 //the number of recovered fish
 init_int n_recovered
 !!cout << "n_recovered: " << n_recovered << endl;
 //the number of time intervals  with homogenous rates per recovered fish
 init_ivector n_intervals_per_fish(1,n_recovered)
 !!cout << "n_intervals_per_fish(n_recovered): " << n_intervals_per_fish(n_recovered) << endl;
 //the total number of time intervals with homogenous rates for all recovered fish
 init_int recovered_total_number_of_intervals
 //the duration of each interval for all recovered fish
 init_vector recovered_intervals(1,recovered_total_number_of_intervals)
 //the reward type for each recovered fish
 init_ivector recovered_reward_type(1,n_recovered)
 //the calendar period for each interval for all recovered fish
 init_ivector recovered_calendar_period(1,recovered_total_number_of_intervals)
 //the post-tagging nonmixing period for each interval for all recovered fish
 init_ivector recovered_nonmixing_period(1,recovered_total_number_of_intervals)
 //the post-tagging tag-induced mortality period for each interval for all recovered fish
 init_ivector recovered_tag_mort_period(1,recovered_total_number_of_intervals)
 //the initial state for each recovered fish
 init_ivector recovered_initial_state(1,n_recovered)
 //the state just prior to capture for each recovered fish
 init_ivector recovered_last_living_state(1,n_recovered)
 //the recovery states for each recovered fish
// init_ivector recovered_end_state(1,n_recovered)
// !!cout << "recovered_last_living_state(n_recovered): " << recovered_last_living_state(n_recovered) << endl;
 
 
 //these determine the number of each type of rate parameter to estimate
 init_int n_mus
 init_int n_Ms
 init_int n_Fs
 init_int n_Ls
 init_int n_rhos
 init_int n_gammas
 init_int n_taus
 !!cout << "n_taus: " << n_taus << endl;

 //these determine the estimation phase of rate parameter
 init_int mu_phase
 init_int M_phase
 init_int F_phase
 init_int L_phase
 init_int rho_phase
 init_int gamma_phase
 init_int tau_phase

 //these determine the number of each type of rate to populate a given instantaneous rate matrix
 init_int n_mus_amat
 //this matrix of indices (1 to n_mus_amat) tells where the mus go in A
 init_imatrix mus_locate(1,n_R,1,n_R)
 
 init_int n_Ms_amat //should be n_R
 init_int n_Fs_amat //should be n_R
 init_int n_Ls_amat //should be n_R
 init_int n_rhos_amat //should be n_R
 init_int n_gammas_amat //should be n_R
 init_int n_taus_amat //should be n_R

 !!cout << "n_taus_amat: " << n_taus_amat << endl;

 //indicates which estimated mu parameter is associated with mus needed for period-specific instantaneous rate matrices
 init_imatrix mus_index(1,n_calendar_periods,1,n_mus_amat)
  !!cout << "mus_index: " << endl;
  !!cout << mus_index << endl;
 //indicates which estimated M parameter is associated with Ms needed for period-specific instantaneous rate matrices
 init_imatrix Ms_index(1,n_calendar_periods,1,n_Ms_amat)
  !!cout << "Ms_index: " << endl;
  !!cout << Ms_index << endl;
 //indicates which estimated F parameter is associated with Fs needed for period-specific instantaneous rate matrices
 init_imatrix Fs_index(1,n_calendar_periods,1,n_Fs_amat)
  !!cout << "Fs_index: " << endl;
  !!cout << Fs_index << endl;
 //indicates which estimated L parameter (tag shedding) is associated with Ls needed for period-specific instantaneous rate matrices
 init_imatrix Ls_index(1,n_calendar_periods,1,n_Ls_amat)
  !!cout << "Ls_index: " << endl;
  !!cout << Ls_index << endl;
 //indicates which estimated rho parameter (reporting prob) is associated with rhos needed for specific reward type
 init_imatrix rhos_index(1,n_reward_types,1,n_rhos_amat)
  !!cout << "rhos_index: " << endl;
  !!cout << rhos_index << endl;
 //indicates which estimated gamma parameter (nonmixing scalar) is associated with gammas needed for specific non-mixing period
 init_imatrix gammas_index(1,n_nonmixing_periods,1,n_gammas_amat)
  !!cout << "gammas_index: " << endl;
  !!cout << gammas_index << endl;
 //indicates which estimated tau parameter (tag induced mort) is associated with taus needed for specific tag-induced mort period
 init_imatrix taus_index(1,n_tag_mort_periods,1,n_taus_amat)

  !!cout << "taus_index: " << endl;
  !!cout << taus_index << endl;
 
 //number of states (2*n_R + 2)
 int n_S
 //number of states where fish is alive (2*n_R)
 int n_S_alive
 //number of states where fish is not reported (2*n_R + 1)
 int n_S_nonreport
 //the states corresponding to being alive (1,2,...,2*n_R)
 ivector living_states(1,n_S_alive) 
 //the states corresponding to being caught and reported (2*n_R+1)
 int report_state
 //the states corresponding to not being not reported (1,2,...,2*n_R, 2*n_R+2)
 ivector non_report_states(1,2*n_R+1)
 
 LOCAL_CALCS
  n_S_alive = 2 * n_R;
  cout << "n_S_alive = 2 * n_R = " << n_S_alive << endl;
  n_S = n_S_alive + 2;
  cout << "n_S = 2 * n_R + 2 = " << n_S << endl;
  report_state = n_S_alive + 1;
  cout << "report_state = n_S_alive + 2 = " << report_state << endl;
  for(ipredict = 1; ipredict <= n_S_alive; ipredict++) non_report_states(ipredict) = ipredict;
  non_report_states(n_S_alive + 1) = n_S_alive + 2;
  cout << "non_report_states: " << non_report_states << endl;
 END_CALCS
 
INITIALIZATION_SECTION
 //PROVIDE INITIAL PARAMETER VALUES
 //best to just supply a .pin file
 log_mus -2.302585
 log_Fs -1.609438
 log_Ms -1.609438
 log_Ls -2.302585
 log_gammas 0.0
 log_taus -3.0
 logit_rhos -1.386294
 
PARAMETER_SECTION
 //DECLARE MODEL PARAMETERS AND VARIABLES

 init_vector log_mus(1,n_mus,mu_phase)
 init_vector log_Ms(1,n_Ms,M_phase)
 init_vector log_Fs(1,n_Fs,F_phase)
 init_vector log_Ls(1,n_Ls,L_phase)
 init_vector logit_rhos(1,n_rhos,rho_phase)
 init_vector log_gammas(1,n_gammas,gamma_phase)
 init_vector log_taus(1,n_taus,tau_phase)
 
 vector mus(1,n_mus)
 vector Ms(1,n_Ms)
 vector Fs(1,n_Fs)
 vector Ls(1,n_Ls)
 vector rhos(1,n_rhos)
 vector gammas(1,n_gammas)
 vector taus(1,n_taus)
  
 vector ll_recovered(1,n_recovered)
 matrix ll_unrecovered(1,n_release_groups_unrecovered,1,2*n_R)
 3darray predicted_ptm(1,n_release_groups_unrecovered,1,3*n_R+1,1,3*n_R+1)
 
 objective_function_value obj_val

RUNTIME_SECTION
 //no tabs!!
 convergence_criteria 1e-4;

PRELIMINARY_CALCS_SECTION

PROCEDURE_SECTION

	obj_val = 0.0;
	ll_unrecovered = 0.0;
	ll_recovered = 0.0;
	
	mus = mfexp(log_mus);
	Fs = mfexp(log_Fs);
	Ms = mfexp(log_Ms);
	Ls = mfexp(log_Ls);
	gammas = mfexp(log_gammas);
	rhos = 1/(1+mfexp(-logit_rhos));
	taus = mfexp(log_taus);	
	
//	cout << "mus: " << mus << endl;
//	cout << "Ms: " << Ms << endl;
//	cout << "Fs: " << Fs << endl;
//	cout << "Ls: " << Ls << endl;
//	cout << "rhos: " << rhos << endl;
//	cout << "gammas: " << gammas << endl;
//	cout << "taus: " << taus << endl;
//	ad_exit(1);
	
	get_ll_fish_unrecovered();
	get_ll_fish_recovered();
	
	obj_val = -sum(ll_recovered) - sum(ll_unrecovered);
//	cout << "obj_val: " << obj_val << endl;
//	ad_exit(1);

FUNCTION get_ll_fish_unrecovered
	
	dvar_matrix P(1,n_S,1,n_S);
	int k = 1;	
//	cout <<"here" << endl;
	dvar_vector temp_ps(1,n_S);
	for(int i = 1; i <= n_release_groups_unrecovered; i++){
		P = identity_matrix(1,n_S);
		for(int j = 1; j <= n_intervals_per_group(i); j++){
			P = P * expm(get_A(unrecovered_calendar_period(k), unrecovered_reward_type(i), unrecovered_nonmixing_period(k), 
				unrecovered_tag_mort_period(k)) * unrecovered_intervals(k));
			k += 1;
		}
//		cout <<"here" << endl;
		for(int l = 1; l <= n_S_alive; l++){
			if(unrecovered_n_released_by_state(i,l) > 0) {
				temp_ps = P(l);
				ll_unrecovered(i,l) =  unrecovered_n_released_by_state(i,l) * log(sum(temp_ps(non_report_states)));
			}
		}
//	cout <<"here" << endl;
	}
//	cout <<"here" << endl;
	if(last_phase()){
		kpredict = 1;	
		for(ipredict = 1; ipredict <= n_release_groups_unrecovered; ipredict++){
//		cout <<"ipredict: "<< ipredict << endl;
			predicted_ptm(ipredict) = identity_matrix(1,3*n_R+1);
			for(int jpredict = 1; jpredict <= n_intervals_per_group(ipredict); jpredict++){
//		cout <<"jpredict: "<< jpredict << endl;
				predicted_ptm(ipredict) = predicted_ptm(ipredict) * 
					expm(get_A_larger(unrecovered_calendar_period(kpredict), unrecovered_reward_type(ipredict), 
					unrecovered_nonmixing_period(kpredict), unrecovered_tag_mort_period(kpredict)) * 
					unrecovered_intervals(kpredict));
				kpredict += 1;	
			}
//		cout << "predicted_ptm(1): " << endl;
//		cout << predicted_ptm(1) << endl;
//		ad_exit(1);		
		}
	}
//	cout << "ll_unrecovered: " << endl;
//	cout << ll_unrecovered << endl;
//	ad_exit(1);

FUNCTION get_ll_fish_recovered
	
	dvar_matrix P(1,n_S,1,n_S);
	int k = 1;	
	
	for(int i = 1; i <= n_recovered; i++){
		P = identity_matrix(1,n_S);
		for(int j = 1; j <= n_intervals_per_fish(i); j++){
			P = P * expm(get_A(recovered_calendar_period(k),recovered_reward_type(i), recovered_nonmixing_period(k), 
				recovered_tag_mort_period(k)) * recovered_intervals(k));
			k += 1;
		}
		ll_recovered(i) = log(P(recovered_initial_state(i), recovered_last_living_state(i))) + 
			log(get_A(recovered_calendar_period(k-1),recovered_reward_type(i), recovered_nonmixing_period(k-1), 
				recovered_tag_mort_period(k-1))(recovered_last_living_state(i), report_state));
//*		if(i ==3){
//			cout << "P:" << endl;
//			cout << P << endl;
//			cout << "A: " << get_A(recovered_calendar_period(k-1),recovered_reward_type(i), recovered_nonmixing_period(k-1), 
//				recovered_tag_mort_period(k-1)) << endl;
//			cout << "reward_type: " << recovered_reward_type(i) << endl;
//			cout << "nonmixing_period: " << recovered_nonmixing_period(k-1) << endl;
//			cout << "tag_mort_period: " << recovered_tag_mort_period(k-1) << endl;
//			cout << "initial_state: " << recovered_initial_state(i) << endl;
//			cout << "last_living_state: " << recovered_last_living_state(i) << endl;
//			cout << "P(initial_state,last_living_state): " << P(recovered_initial_state(i), recovered_last_living_state(i)) << endl;
//			cout << "final a_ij: " << get_A(recovered_calendar_period(k-1),recovered_reward_type(i), recovered_nonmixing_period(k-1), 
//				recovered_tag_mort_period(k-1))(recovered_last_living_state(i), report_state) << endl;
//			cout << "ll_recovered(3): " << ll_recovered(3) << endl;
//			
//			ad_exit(1);
//		}

	}
//	cout << "ll_recovered(1): " << ll_recovered(1) << endl;
//	cout << "sum(ll_recovered): " << sum(ll_recovered) << endl;
		

FUNCTION dvar_matrix get_A(int& calendar_period, int& reward_type, int& nonmixing_period, int& tag_mort_period)
	dvar_vector temp_rhos(1,n_rhos_amat);
	dvar_vector temp_gammas(1,n_gammas_amat);
	dvar_vector temp_taus(1,n_taus_amat);
	dvar_matrix A(1,n_S,1,n_S);
	
//	cout << "rhos_index: " << endl;
//	cout << rhos_index << endl;
//	cout << "reward_type: " << endl;
//	cout << reward_type << endl;
	if(reward_type == 0) temp_rhos = 1.0; //highest reward type
	else temp_rhos = rhos(rhos_index(reward_type)); 
//	cout << "temp_rhos: " << endl;
//	cout << temp_rhos << endl;
	
	if(nonmixing_period == 0) temp_gammas = 1.0; //fully mixed
	else temp_gammas = gammas(gammas_index(nonmixing_period));
	if(tag_mort_period == 0) temp_taus = 0.0; //no more tag-induced mortality
	else temp_taus = taus(taus_index(tag_mort_period));
	//A = make_A(mus(mus_index(calendar_period)), Fs(Fs_index(calendar_period)), Ms(Ms_index(calendar_period)), Ls(Ls_index(calendar_period)), 
	//	temp_rhos, temp_gammas, temp_taus, n_R);
	A = make_A_smaller(mus(mus_index(calendar_period)), Fs(Fs_index(calendar_period)), Ms(Ms_index(calendar_period)), Ls(Ls_index(calendar_period)), 
		temp_rhos, temp_gammas, temp_taus, n_R);
	return(A);

FUNCTION dvar_matrix get_A_larger(int& calendar_period, int& reward_type, int& nonmixing_period, int& tag_mort_period)
	dvar_vector temp_rhos(1,n_rhos_amat);
	dvar_vector temp_gammas(1,n_gammas_amat);
	dvar_vector temp_taus(1,n_taus_amat);
	dvar_matrix A(1,3*n_R+1,1,3*n_R+1);
	
	if(reward_type == 0) temp_rhos = 1.0; //highest reward type
	else temp_rhos = rhos(rhos_index(reward_type)); 
	
	if(nonmixing_period == 0) temp_gammas = 1.0; //fully mixed
	else temp_gammas = gammas(gammas_index(nonmixing_period));
	if(tag_mort_period == 0) temp_taus = 0.0; //no more tag-induced mortality
	else temp_taus = taus(taus_index(tag_mort_period));
	A = make_A_larger(mus(mus_index(calendar_period)), Fs(Fs_index(calendar_period)), Ms(Ms_index(calendar_period)), Ls(Ls_index(calendar_period)), 
		temp_rhos, temp_gammas, temp_taus, n_R);
	return(A);

FUNCTION dvar_matrix make_A(dvar_vector mus, dvar_vector Fs, dvar_vector Ms, dvar_vector Ls, dvar_vector rhos, dvar_vector gammas, dvar_vector taus, int n_R)
	
	//this will incorporate tag shedding, tag-induced mortality, reporting rate, and non-mixing for 3 regions that are all connected
	
	//states 1,...,n_R correspond to being alive with 2 tags in a given region
	//states n_R + 1,...,2*n_R correspond to being alive with 1 tag in a given region
	//state 2*n_R + 1 corresponds to being alive with no tags
	//state 2*n_R + 2 corresponds to caught and reported
	//state 2*n_R + 3 corresponds to caught and not reported
	//state 2*n_R + 4 corresponds to death from other causes
	
	//mu_ij migration rates
	//F_1, ..., F_n_R, fishing mortality 
	//M_1, ..., M_n_R, natural mortality
	//L_1, ..., L_n_R, tag-shedding
	//rho_1, ..., rho_n_R, reporting probability
	//gamma_1, ..., gamma_n_R, non-mixing scalar
	//tau_1, ..., tau_n_R, tag-induced mortality
	
	//the K x K instantaneous rate matrix (A) returned is (K = 2*n_R+4, n_R is number of regions)
	//    	          1, ...,    n_R,   n_R + 1, ...,     2*n_R, 2*n_R + 1,               2*n_R + 2,                   2*n_R + 3,      2*n_R + 4
	//1		        a_1, ..., m_1n_R,     2*L_1, ...,         0,         0,       rep_1*gamma_1*F_1,        (1-rep_1)*gamma1*F_1,    tau_1 + M_1
	//...           ..., ...,    ...,       ..., ...,       ...,       ...,                     ...,                         ...,            ...
	//n_R       mu_n_R1, ...,  a_n_R,         0, ...,   2*L_n_R,         0, rep_n_R*gamma_n_R*F_n_R, (1-rep_n_R)*gamma_n_R*F_n_R, tau_n_R +M_n_R
	//n_R + 1         0, ...,      0, a_(n_R+1), ...,   mu_1n_R,       L_1,       rep_1*gamma_1*F_1,        (1-rep_1)*gamma1*F_1,    tau_1 + M_1
	//...           ..., ...,    ...,       ..., ...,       ...,       ...,                     ...,                         ...,            ...
	//2*n_R           0, ...,      0,   mu_n_R1, ..., a_(2*n_R),     L_n_R, rep_n_R*gamma_n_R*F_n_R, (1-rep_n_R)*gamma_n_R*F_n_R, tau_n_R +M_n_R
	//2*n_R + 1       0, ...,      0,         0, ...,         0,         0,                       0,                           0,              0
	//...           ..., ...,    ...,       ..., ...,       ...,       ...,                     ...,                         ...,            ...
	//2*n_R + 4       0, ...,      0,         0, ...,         0,         0,                       0,                           0,              0
	// all zeros 4 rows below 2*n_R!!!!
	
	dvar_matrix A(1,n_S,1,n_S);
	A = identity_matrix(1,n_S);
		
	for(int i = 1; i <= n_R; i++){
		//reported and caught
		A(i,2*n_R + 2) = A(i+n_R,2*n_R + 2) = rhos(i)*gammas(i)*Fs(i);
		//unreported and caught
		A(i,2*n_R + 3) = A(i+n_R,2*n_R + 3) = (1-rhos(i))*gammas(i)*Fs(i);
		//natural mortality
		A(i,2*n_R + 4) = A(i+n_R,2*n_R + 4) = Ms(i) + taus(i);
		//starts with two tags and loses one
		A(i,i+n_R) = 2*Ls(i);
		//starts with one tag and loses it
		A(i+n_R,2*n_R + 1) = Ls(i);
		for(int j = 1; j <= n_R; j++){ //put the migration rate parameters in the right positions of A
			if(mus_locate(i,j)>0) A(i,j) = A(i+n_R,j+n_R)= mus(mus_locate(i,j));
		}
		//need to make the diagonal equal to -1*sum of the other elements in the row for the first 2*n_R rows
		A(i,i) = -(sum(A(i)) - A(i,i)); 
		A(i+n_R,i+n_R) = -(sum(A(i+n_R)) - A(i+n_R,i+n_R)); 
	}
	for(int i = 2*n_R + 1; i <= 2*n_R + 4; i++) A(i,i) = 0.0; //fill in the rest of the diagonals with 0
		
	return A;


FUNCTION dvar_matrix make_A_smaller(dvar_vector mus, dvar_vector Fs, dvar_vector Ms, dvar_vector Ls, dvar_vector rhos, dvar_vector gammas, dvar_vector taus, int n_R)
	
	//this will incorporate tag shedding, tag-induced mortality, reporting rate, and non-mixing for 3 regions that are all connected
	
	//states 1,...,n_R correspond to being alive with 2 tags in a given region
	//states n_R + 1,...,2*n_R correspond to being alive with 1 tag in a given region
	//state 2*n_R + 1 corresponds to caught and reported
	//state 2*n_R + 2 corresponds to being alive with no tags, caught and not reported, or death from other causes
	
	//mu_ij migration rates
	//F_1, ..., F_n_R, fishing mortality 
	//M_1, ..., M_n_R, natural mortality
	//L_1, ..., L_n_R, tag-shedding
	//rho_1, ..., rho_n_R, reporting probability
	//gamma_1, ..., gamma_n_R, non-mixing scalar
	//tau_1, ..., tau_n_R, tag-induced mortality
	
	//the K x K instantaneous rate matrix (A) returned is (K = 2*n_R+2, n_R is number of regions)
	//    	          1, ...,    n_R,   n_R + 1, ...,     2*n_R,               2*n_R + 2, 	    	    		   				      2*n_R + 2
	//1		        a_1, ..., m_1n_R,     2*L_1, ...,         0,       rep_1*gamma_1*F_1,                0 + (1-rep_1)*gamma1*F_1 + tau_1 + M_1
	//...           ..., ...,    ...,       ..., ...,       ...,                     ...,         								            ...
	//n_R       mu_n_R1, ...,  a_n_R,         0, ...,   2*L_n_R, rep_n_R*gamma_n_R*F_n_R,     0 + (1-rep_n_R)*gamma_n_R*F_n_R + tau_n_R + M_n_R
	//n_R + 1         0, ...,      0, a_(n_R+1), ...,   mu_1n_R,       rep_1*gamma_1*F_1,              L_1 + (1-rep_1)*gamma1*F_1 + tau_1 + M_1
	//...           ..., ...,    ...,       ..., ...,       ...,                     ...,                                                   ...
	//2*n_R           0, ...,      0,   mu_n_R1, ..., a_(2*n_R), rep_n_R*gamma_n_R*F_n_R, L_n_R + (1-rep_n_R)*gamma_n_R*F_n_R + tau_n_R + M_n_R
	//2*n_R + 1       0, ...,      0,         0, ...,         0,                       0,         0,                           0,              0
	//...           ..., ...,    ...,       ..., ...,       ...,                     ...,       ...,                         ...,            ...
	//2*n_R + 4       0, ...,      0,         0, ...,         0,                       0,         0,                           0,              0
	// all zeros 4 rows below 2*n_R!!!!
	
	dvar_matrix A(1,2*n_R + 2,1,2*n_R + 2);
	A = identity_matrix(1,2*n_R + 2);
//	cout << "make_A_smaller" << endl;
//	cout << "mus: " << mus << endl;
//	cout << "Ms: " << Ms << endl;
//	cout << "Fs: " << Fs << endl;
//	cout << "Ls: " << Ls << endl;
//	cout << "rhos: " << rhos << endl;
//	cout << "gammas: " << gammas << endl;
//	cout << "taus: " << taus << endl;
		
	for(int i = 1; i <= n_R; i++){
		//reported and caught
		A(i,2*n_R + 1) = A(i+n_R,2*n_R + 1) = rhos(i)*gammas(i)*Fs(i);
		//unreported and caught + natural mortality + tag-induced mortality
		A(i,2*n_R + 2) = A(i+n_R,2*n_R + 2) = (1-rhos(i))*gammas(i)*Fs(i) + Ms(i) + taus(i);
		//starts with two tags and loses one
		A(i,i+n_R) = 2*Ls(i);
		//starts with one tag and loses it (add it to unreported and caught + natural mortality + tag-induced mortality)
		A(i+n_R,2*n_R + 2) += Ls(i);
		for(int j = 1; j <= n_R; j++){ //put the migration rate parameters in the right positions of A
			if(mus_locate(i,j)>0) A(i,j) = A(i+n_R,j+n_R)= mus(mus_locate(i,j));
		}
		//need to make the diagonal equal to -1*sum of the other elements in the row for the first 2*n_R rows
		A(i,i) = -(sum(A(i)) - A(i,i)); 
		A(i+n_R,i+n_R) = -(sum(A(i+n_R)) - A(i+n_R,i+n_R)); 
	}
	for(int i = 2*n_R + 1; i <= 2*n_R + 2; i++) A(i,i) = 0.0; //fill in the rest of the diagonals with 0
		
	return A;

FUNCTION dvar_matrix make_A_larger(dvar_vector mus, dvar_vector Fs, dvar_vector Ms, dvar_vector Ls, dvar_vector rhos, dvar_vector gammas, dvar_vector taus, int n_R)
	
	//this will incorporate tag shedding, tag-induced mortality, reporting rate, and non-mixing for 3 regions that are all connected
	
	//states 1,...,n_R correspond to being alive with 2 tags in a given region
	//states n_R + 1,...,2*n_R correspond to being alive with 1 tag in a given region
	//state 2*n_R + 1 corresponds to caught and reported in region 1
	//state 2*n_R + ... corresponds to caught and reported in region ...
	//state 2*n_R + n_R corresponds to caught and reported in region n_R
	//state 3*n_R + 1 corresponds to being alive with no tags, caught and not reported, or death from other causes
		
	dvar_matrix A(1,3*n_R+1,1,3*n_R+1);
	A = identity_matrix(1,3*n_R+1);
//	cout << "mus: " << mus << endl;
//	cout << "Ms: " << Ms << endl;
//	cout << "Fs: " << Fs << endl;
//	cout << "Ls: " << Ls << endl;
//	cout << "rhos: " << rhos << endl;
//	cout << "gammas: " << gammas << endl;
//	cout << "taus: " << taus << endl;
		
	for(int i = 1; i <= n_R; i++){
		//starts with two tags and loses one
		A(i,i+n_R) = 2*Ls(i);
		//reported and caught
		A(i,2*n_R + i) = A(i+n_R,2*n_R + i) = rhos(i)*gammas(i)*Fs(i);
		//starts with two tags and gets caught and not reported, or dies due to other causes (nat and tagging mortality)
		A(i,3*n_R + 1) = (1-rhos(i))*gammas(i)*Fs(i) + Ms(i) + taus(i);
		//starts with one tag and loses it, gets caught and not reported, or dies due to other causes (nat and tagging mortality)
		A(i+n_R,3*n_R + 1) = Ls(i) + (1-rhos(i))*gammas(i)*Fs(i) + Ms(i) + taus(i);
		for(int j = 1; j <= n_R; j++){ //put the migration rate parameters in the right positions of A
			if(mus_locate(i,j)>0) A(i,j) = A(i+n_R,j+n_R)= mus(mus_locate(i,j));
		}
		//need to make the diagonal equal to -1*sum of the other elements in the row for the first 2*n_R rows
		A(i,i) = -(sum(A(i)) - A(i,i)); 
		A(i+n_R,i+n_R) = -(sum(A(i+n_R)) - A(i+n_R,i+n_R)); 
	}
	for(int i = 2*n_R + 1; i <= 3*n_R + 1; i++) A(i,i) = 0.0; //fill in the rest of the diagonals with 0
//	cout << "A: " << endl;
//	cout << A << endl;
	return A;



REPORT_SECTION
  //OUTPUT RESULTS TO FILE "fsct_tr_1.rep"

  report << "finite-state continuous-time tagging model" << endl;
  report << "OBJECTIVE FUNCTION VALUE: " << endl;
  report << obj_val << endl;

  report << "pars:" << endl;
  report << log_mus << " " << log_Ms << " " << log_Fs << " " << log_Ls << " " << logit_rhos << " " << log_gammas 
  	<< " " << log_taus << endl;

  report << "ll_recovered: " << endl;
  report << sum(ll_recovered) << endl;

  report << "ll_unrecovered: " << endl;
  report << sum(ll_unrecovered) << endl;

  report << "mus:" << endl;
  report << mus << endl;

  report << "Ms:" << endl;
  report << Ms << endl;
  
  report << "Fs:" << endl;
  report << Fs << endl;

  report << "Ls:" << endl;
  report << Ls << endl;
  
  report << "rhos:" << endl;
  report << rhos << endl;

  report << "gammas:" << endl;
  report << gammas << endl;

  report << "taus:" << endl;
  report << taus << endl;
  
  for(ipredict = 1; ipredict <= n_release_groups_unrecovered; ipredict++){
    report << "predicted_ptm(" << ipredict << "):" << endl;
	report << predicted_ptm(ipredict) << endl;
  }

  report << "ll_recovered: " << endl;
  report << ll_recovered << endl;

FINAL_SECTION
  //Calculates how long is taking to run
  // this code is based on the Widow Rockfish model (from Erik H. Williams, NMFS-Santa Cruz, now Beaufort)  
  time(&finish);
  elapsed_time = difftime(finish,start);
  hour = long(elapsed_time)/3600;
  minute = long(elapsed_time)%3600/60;
  second = (long(elapsed_time)%3600)%60;
  cout<<endl<<endl<<"starting time: "<<ctime(&start);
  cout<<"finishing time: "<<ctime(&finish);
  cout<<"This run took: ";
  cout<<hour<<" hours, "<<minute<<" minutes, "<<second<<" seconds."<<endl<<endl<<endl;
